Skip to main content

Full text of "Image Deblurring Using Derivative Compressed Sensing for Optical Imaging Application"

See other formats


Image Deblurring Using Derivative Compressed 
Sensing for Optical Imaging Application 

Mohammad Rostami, Student Member, IEEE, Oleg Michailovich, Member, IEEE, and Zhou Wang, Member, IEEE 










Abstract — Reconstruction of multidimensional signals from the 
samples of their partial derivatives is known to be a standard 
problem in inverse theory. Such and similar problems routinely 
arise in numerous areas of applied sciences, including optical 
imaging, laser interferometry, computer vision, remote sensing 
and control. Though being ill-posed in nature, the above problem 
can be solved in a unique and stable manner, provided proper 
regularization and relevant boundary conditions. In this paper, 
however, a more challenging setup is addressed, in which one 
has to recover an image of interest from its noisy and blurry 
version, while the only information available about the imaging 
system at hand is the amplitude of the generalized pupil function 
(GPF) along with partial observations of the gradient of GPF's 
phase. In this case, the phase-related information is collected 
using a simplified version of the Shack-Hartmann interferometer, 
followed by recovering the entire phase by means of derivative 
compressed sensing. Subsequently, the estimated phase can be 
combined with the amplitude of the GPF to produce an estimate 
of the point spread function (PSF), whose knowledge is essential 
for subsequent image deconvolution. In summary, the principal 
contribution of this work is twofold. First, we demonstrate how to 
simplify the construction of the Shack-Hartmann interferometer 
so as to make it less expensive and hence more accessible. Second, 
it is shown by means of numerical experiments that the above 
simplification and its associated solution scheme produce image 
reconstructions of the quality comparable to those obtained using 
dense sampling of the GPF phase. 

Index Terms — Deconvolution, inverse problem, derivative com- 
pressive sampling, and Shack-Hartmann interferometer. 

I. Introduction 

THE necessity to recover digital images from their dis- 
torted and noisy observations arises in a multitude of 
practical applications, with some specific examples including 
image denoising, super-resolution, image restoration, and wa- 
termarking, just to name a few |[l)-|j4J. In such cases, it is 
standard to assume that the observed image v is formed by 
subjecting the original image u to convolution with a point 
spread functiorT] (PSF) i, followed by contamination by white 
Gaussian noise (WGN) v. Thus, formally. 

i-^ u + V. 


While u and v can be regarded as general members of the 
signal space L2(ri) of real-valued functions on 51 C M?, the 
PSF i is normally a much smoother function, with effectively 
band-limited spectrum. As a result, the convolution with i has 
a destructive effect on the informational content of u, in which 

All the authors are with the Department of Electrical and Computer 
Engineering, University of Waterloo, N2L 3G1 Ontario. 

' Note that, in optical imaging, this function is also referred to as an impulse 
transfer function |5 1. 

case V typically has a substantially reduced set of features with 
respect to u. This makes the problem of reconstruction of u 
from V a problem of significant practical importance 161. 

Reconstruction of the original image u from v can be carried 
out within the framework of image deconvolution, which is a 
specific instance of a more general class of inverse problems 
JT). Most of such methods are Bayesian in nature, in which 
case the information lost in the process of convolution with i 
is recovered by requiring the optimal solution to reside within 
a predefined functional class |8 1, |9|. Thus, for example, in the 
case when u is known to be an image of bounded variation, the 
above regularization leads to the famous Rudin-Osher-Fatemi 
reconstruction scheme, in which u is estimated as a solution 

to the following problem 1 10 1, 1 1 1 1 

arg mm 



|Vu| dxdy > 


where a > is the regularization parameter. It should be noted 
that, if the PSF obeys J i dxdy ^ 0, the problem O) is strictly 
convex and therefore admits a unique minimizer, which can 
be computed by a spectrum of available algorithms pO) , pTj . 

A particularly non-trivial version of deconvolution is com- 
monly referred to as blind. In this case, the original image 
u is to be estimated without the knowledge of the PSF I?). 
In this paper, however, we follow the philosophy of hybrid 
deconvolution |12|, which takes advantage of any partial 
information on the PSF to improve the image reconstruction. 
Thus, in the algorithm described in this paper, the original 
image u will be recovered from v and some partial information 
on i. 

Optical (and, in particular, turbulent) imaging is unarguably 
the field of applied sciences from which the notion of decon- 
volution has originally emanated p3j-p7j. In short-exposure 
imaging, however, computational methods of image restoration 
are still superseded by adaptive optics. As recently as a 
decade ago, the use of adaptive optics would have been 
considered as the only practical option. Nowadays, however, 
with the advent of distributed cluster computing and GPU- 
based image processing, it seems to be time to revisit the cost- 
to-performance characteristics of the existing tools of adaptive 
optics. Thus, in this work, our focus is on a specific tool of 
adaptive optics, known as the Shack-Hartmann interferometer 
1 16 1, p7| . Instead of completely excluding the interferometer 
from our measurement system, we propose to modify its 
construction through reducing the number of its local wave- 
front lenses. Although the advantages of such a simplification 
are immediate to see, its main shortcoming is obvious as 
well: the smaller the number of lenses is, the stronger is 


the effect of undersampling and aliasing. Accordingly, to 
overcome this problem, we propose to augment the modified 
Shack-Hartmann interferometer by subjecting its output to the 

derivative compressed sensing (DCS) algorithm of 1 18 1. As it 
will be shown later in the paper, the PSF i is determined by 
a generalized pupil function P, which can be expressed in a 
polar form as P = A e^"^. While the amplitude A can be mea- 
sured via calibration or computed as a function of the aperture 
geometry, the phase is often influenced by environmental 
effects and hence it needs to be recovered from observations. 
It will be shown below that DCS is particularly well suited 
for reconstruction of </) from incomplete measurements of 
its partial differences. Such an estimate can be subsequently 
combined with A to yield an estimate of the PSF i, which 
can in turn be used by a deconvolution algorithm. Thus, the 
proposed method for estimation of the PSF and subsequent 
deconvolution of u can be regarded as a hybrid deconvolution 
technique, which comes to simplify the design and complexity 
of adaptive optics on the one hand, and to make the process 
of reconstruction of optical images as automatic as possible, 
on the other hand. 

The rest of the paper is organized as follows. Section II 
summarizes basic technical preliminaries. In Section III, we 
describe the SH interferometer as well as phase measurements 
in optical imaging. In Section IV, we explain DCS and our new 
approach to solve it. In Section V we describe deconvolution 
process to recover the original image. Experimental results are 
presented in Section VI, while Section VII finahzes the paper 
with a discussion and main conclusions. 

II. Technical Preliminaries 

In short exposure imaging, due to aberrations in the imaging 
system induced by, e.g., atmospheric turbulence, the impulse 
response of an optical imaging system is often unknown p9) . 
In order to better understand the setup under consideration, we 
first note that, in optical imaging, the PSF i is obtained from 
an amplitude spread function (ASF) h as i :— |ft,p. The ASF, 
in turn, is defined in terms of the generalized pupil function 
(GPF) P(x, y) as given by 1,20} 






dxdy, (3) 

where Zi is the focal distance and A^, is the optical wavelength. 
Being a complex-valued quantity, P{x, y) can be represented 
in terms of its amplitude A{x, y) and phase (/)(x, y) as 



Here, the GPF amplitude A{x,y) (which is sometimes simply 
referred to as the aperture function) is normally a function 
of the aperture geometry. Thus, for instance, in the case of a 
circular aperture, A{x, y) can be defined as ||19| 

A{r) = 



where D denotes the pupil diameter. Thus, given (j){x,y), one 
could determine h and therefore i. Unfortunately, the phase 
(j){x, y) does not have an analytic expression, and it has to be 

measured in practice using such tools as the Shack-Hartmann 
interferometer (SHI) | fT6l . 

As will be discussed later in the paper, the SHI is capable 
of sensing the partial derivatives of (t>{x,y). Needless to say, 
in order to minimize the effect of aliasing on the estimation 
result, an accurate reconstruction of 0(a:, y) requires taking a 
fairly large number of the samples of V0(x, y) |21 1. In some 
applications, the number of sampling points (as defined by the 
number of local wavefront lenses) reaches as many as a few 
thousands. It goes without saying that reducing the number of 
lenses would have a positive impact on the SHI in terms of its 
cost and approachability. Alas, such a reduction is impossible 
without undersampling, which tends to have formidable effect 
on the overall quality of phase estimation. 

In this paper, to minimize the effect of undersampling, we 
exploit DCS |18|. As opposed to the classical compressed 
sensing (CCS) |22| , in addition to the sparsifing constraints, 
DCS also uses constraints which are intrinsic in the definition 
of partial derivatives. Using these additional constraints - 
which are called the cross-derivative constraints - allows sub- 
stantially improving the quality of reconstruction of (t)[x,y), 
as compared to the case of CCS-based estimation. 

III. Shack-Hartmann Interferometer (SHI) 

As it was mentioned earlier, the SHI is typically used to 
measure the gradient V0(x, y) of the GPF phase 0(a;, y), from 
which the values of the latter can be subsequently estimated. 
To this end, the unknown phase (f){x, y) is assumed to be 
expandable in terms of some basis functions {Zk}'^^^, viz. 

^ ' oo 

(l^{x,y)^^akZk{x,y), (6) 


where the representation coefficients {ak}kLo are assumed 
to be unique and stably computable. Note that, in this case, 
the datum of {flfej^Q uniquely identifies (f){x,y), while the 
coefficients {a/d^g can be estimated due to the linearity of 
(|6]l which suggests 

V(?!)(x, y) = ^ak ^Zk{x, y), 



The most frequent choice of {^fe}^o ^'^ adaptive optics is 
the Zernike polynomials (aka Zemike functions) [20 1. These 
polynomials constitute an orthonormal basis in the space of 
square-integrable functions defined over the unit disk in M^. 
Zernike polynomials can be subdivided in two subsets of the 
even Z™ and odd Z,7™ Zernike polynomials which have very 
convenient analytical definitions as given by 



where m and n are nonnegative integers with n > m, < 
If < 271 is the azimuthal angle, and < p < 1 is the radial 
distance. The radial polynomials _R™ are defined as 



(-!)'=(« -fc)! 

^ A:!((n + m)/2-/s)!((n-m)/2-fc)! 

11-2 k 



Fig. 1. An example of a 10 X 10 SHI array on a circular aperture. The 
shading indicates those blocks (i.e., lenses) which are rendered idle. 

Note that, since the Zernike polynomials above are defined 
using polar coordinates, it makes sense to re-express the 
phase (p and its gradient in the polar coordinate system as 
well. (Technically, this would amount to replacing x and y 
in (|6]l-(j7]i by p and ip, respectively.) Moreover, due to the 
property of the Zernike polynomials to be an orthonormal 
basis, the representation coefficients {a^l^Q in in ((6])-(|7]i can 
be computed by orthogonal projection, namely 


(t){p,(p) Zk{p,(p) pdpdip 


In practice, however, (j){p, (p) is unknown and therefore the 
coefficients {flfej^o n^ed to be estimated by other means. 
Thus, in the case of the SHI, the coefficients can be estimated 
from a finite set of discrete measurements of V0(p, (p). 

The main function of the SHI is to acquire discrete measure- 
ments of V0 by means of linearization. The linearization takes 
advantage of subdividing a (circular) aperture into rectangular 
blocks with their sides formed by a uniform rectangular lattice. 
An example of such a subdivision is shown in Fig. [Tlfor the 
case of a 10 X 10 lattice grid. Subsequently, it is assumed 
that the grid is sufficiently fine to approximate a restriction of 
the phase <j) to each of the above blocks by a linear function. 
This results in a piecewise linear approximation of 0, whose 
accuracy improves asymptotically when the lattice size goes 
to infinit>|^ Formally, let n := {{x,y) e M^ | x^ + y^ < jj2y 
be a circular aperture of radius D and 5 = {(x, y) e M^ | 
max{|a;|, \y\} < D} be a square subset of R-^ such that fi C 5. 
Then, for each polar coordinate (p, ip) £ ft and an N x N grid 
of square blocks of size 2D/N x 2D/N, the phase (j) can be 
expressed as 

4>(x,y) ~ ax + by + c^ (12) 

for all (x, y) in a neighbourhood of (p cos (p,p sin ip). The 
approximation in ([T2| suggests that 

V0(x,y)« [a, 6]^ 


^More rigorously, one can show that, as long as </> is uniformly continuous, 
its piecewise linear approximation converges uniformly, as the grid size goes 
to infinity. 

LoisAnij C CD Clip 

Shack-Haattraan Pattern 

Fig. 2. Basic structure of the SHI and a resulting pattern of the focal points. 

where (•)-^ denotes matrix transposition. While c in ( 12 1 can be 

derived from boundary conditions, coefficients a and b should 
be determined via direct measurements. To this end, the SHI 
is endowed with an array of small focusing lenses, which are 
supported over each of the square blocks of the discrete grid. In 
the absence of phase aberrations, the focal points of the lenses 
are spatially identified and registered using a high-resolution 
CCD detector, whose imaging plane is aligned with the plane 
of the focal points. Then, when the wavefront gets distorted 
as a result of, e.g., atmospheric turbulence, the focal points 
are "pushed" towards new spatial positions, which can also be 
pinpointed by the same detector. The resulting displacements 
can therefore be measured and subsequently related to the 
values of V0 at corresponding points. 

To explain how the above procedure can be performed, 
additional notations are in order. Let ft^ denote a finite set 
of spatial coordinates defined according to 





e n 

-D+^(i + ^], t^O,l,...,N~l (14) 

and xl + yj<D^y 

Note that the set fid can be thought of as a set of the 
spatial coordinates of the geometric centres of the wavefront 
lenses, restricted to the domain of aperture O. Under some 
reasonable assumptions, one can then show 1 ,23) that the focal 
displacement A{x,y) — [Ax{x, y), Ay{x,y)Y' measured at 
some {x, y) G fid is related to the value of V(/)(x, y) as given 

V0(a;, y) « -A^{x, y), V(a;, y) G ^d, (15) 

where / is the focal distance of the wavefront lenses. Such a 
measurement setup is depicted in Fig.|2]along with an example 
of measured focal points. 

Now, given a total of M :— #51^ measurements of V0 
over rijj, one can try to recover a useful approximation of (j) 
over the whole Vl in the form of a projection of 4> on the 
linear subspace spanned by all Zernike polynomials up to the 


order L inclusive. In this case, it is possible to estimate the 
representation coefficients {ak\k=a ^^^ solution of 

conditions of compressed sensing far superseded those of the 

subject to appropriate boundary conditions. It is worthwhile 

noting that (16 1 can be rewritten in a vector-matrix form as 



s.t. a h 0, 


where Z is a 2Af x i + 1 matrix composed of the values 
of the partial derivatives of the Zernike polynomials, d is 
a measurement (column) vector of length 2M, and a — 
[oq, fli, . . . , ai]^ is a vector of the representation coefficients 

of (j). The constraint a )^ in (17 1 is supposed to further 

regularize the solution by forcing a to be a member of some 
convex cone as well. Thus, for example, if the mean value of 

(j) can be assumed to be equal to zero, the solution to ( 17 1 can 
be computed as 

a=Z#d, (18) 

where Z* denotes the pseudo-inverse of Z, whose definition 
is unique and stable as long as the row-rank of Z is greater or 
equal to L + 1 (hence suggesting that 2M > L + 1). Having 
estimated a, the phase can be approximated as 

<P{P,'~P) ~ ^akZk(p,(p). 



A higher accuracy of phase estimation requires using higher- 
order Zernike polynomials, which in turn necessitates a pro- 
portional increase in the number of wavefront lenses. More- 
over, as required by the linearization procedure in the SHI, 
the lenses have to be of a relatively small sizes (sometimes, 
on the order of a few microns), which may lead to the use 
of a few thousand lenses per one interferometer. Accordingly, 
to simplify the construction and to reduce the cost of SHIs, 
we propose to substantially reduce the number of wavefront 
lenses, while compensating for the induced information loss 
through the use of DCS, which is detailed next. 

IV. Derivative Compressive Sampling 

A. Classical Compressed Sensing 

Central to signal processing is the Shannon-Nyquist theorem 
p4) , which specifies conditions on which a band-limited 
signal can be stably and uniquely recovered from its discrete 
measurements. However, in around 2005, a different sampling 
theorem was formed which, in some cases, abrogates the 
fundamentals of its predecessor. This new theory, nowadays 
known as compressed sensing (aka compressive sampling), 
asserts that signals, which admit a sparse representation in 
a predefined basis/frame, can be recovered from their discrete 
measurements, whose number is proportional to the ^o-norm 
of the coefficients of the sparse representation. In such a 
case, the sparser the representation of the signal is, the 
smaller can be the number of measurements required for signal 
reconstruction. As a result, cases are numerous in which the 

Shannon-Nyquist sampling |22|, |25|. 

Despite its widespread success in countless applications, the 
theory of compressed sensing is not entirely free of limitations. 
One of such limitations stems from the necessity to use a non- 
linear decoder. Indeed, while in the case of classical sampling, 
the reconstruction of a time-domain signal is implemented 
through linear interpolation (i.e., linear filtering), in the case 
of compressed sensing, the reconstruction involves solution of 
a convex optimization problem. Specifically, let us consider a 
typical setup of classical compressed sensing (CCS), in which 
y £ M'" represents an observed version of x e M", related 
according to 

y = *x, (20) 

where ^ G M™^" is an observation (sampling) matrix with 
n > m. 

The recovery of x from y based on ( [20| is impossible to 
implement in a unique and stable way, unless it is known that 
X is sparse and hence has a relatively low value of ||x||o. In 
such a case, if the sampling matrix ^ satisfies the restricted 
isometry property (RIP) f22|, p5) with respect to a certain 
class of sparse signals to which x is assumed to belong, then 
CCS recovers x as a solution to | |26l , | |27l 

X = argmin{||x'||i | ^x' = y} , (21) 


which is a convex minimization problem, which is straightfor- 
ward to reformulate in terms of linear programming. Moreover, 
in the case when the measurements y are error-prone, a more 
robust version of CCS is to recover x as given by 

X — arg mm 






where e > is a parameter controlling the size of the noise. 
Moreover, it was shown in 1221, ||25l, that the estimation error 

in the signal reconstructed according to ( 22 1 can be bounded 
by a linear function of e. This implies robustness of the CCS 
reconstruction towards the presence of measurement noise. 
It should be finally noted that the optimization problem 

(22 1 can be reformulated in its alternative Lagrangian form. 

in which case one can find 

X — arg mm 





where A > is an optimal Lagrange multiplier p6) . In what 
follows, it is assumed that an optimal value of A is known. 
(For more details on this subject, the reader is referred to p6| 
as well as to the later sections of this paper). 

B. Derivative Compressed Sensing (DCS) 

Let the partial derivatives of evaluated (by means of the 
SHI) at the points of set ilij be column-stacked into vectors 
fj; and iy of length M = i^fld- In what follows, the partial 
derivatives i^ and fj, are assumed to be sparsely representable 
by an orthonormal basis in M*^. Representing such a basis by 
an M X M unitary matrix W, the above assumption suggests 
the existence of two sparse vectors c^, and Cy such that 
Wcx and Wcy provide accurate approximations of the partial 
derivatives of the original phase (j) evaluated over fid- 


Now, the simplification of the SHI proposed in the current 
paper amounts to reducing the number of wavefront lenses 
to a minimum. Formally, such a reduction can be described 
by two n X M sub-sampling matrices '^^ and '^y, where 
n < M. Specifically, let h^ '■— 'i'xix and bj, :— 'i'yiy be 
incomplete (partial) observations of i^ and fj^, respectively. 
Then, the noise-free counterparts of the partial derivatives can 
be approximated by Wc^ and Wc* respectively, where 

< = argmin |i||*,W^< - b,||2 + A,||<||i| (24) 

where p^*) is a vector of Bregman variables (or, equivalently, 
augmented Lagrange multipliers) and (5 > is a user-defined 


argmin<; -W'i'yWc'y - hy\\l 




for some Aa;,Ay > 0. Moreover, in the case when A^^ — Xy, 

the above estimates can be combined together To this end, let 

c = [Cx,Cyf, b = [K,hyf, and A = diag{*^M^, *yTy} e 
M2"x2M^ Then, 

argmm ■^ -||Ac' 

b||2 + A||c'||i 


The c-update step in ( 32 1 has the format of a standard 
basis pursuit de-noising (BPDN) problem f29l, which can 
be solved by a variety of optimization methods [301. In the 
present paper, we used the FISTA algorithm of |31) due to the 
simplicity of its implementation as well as for its remarkable 
convergence properties. It should be noted that the algorithm 
does not require explicitly defining the matrices A and B. 
Only the operations of multiplication by these matrices and 
their transposes need to be known, which can be implemented 
in an implicit and computationally efficient manner 

Once an optimal c* is recovered, it can be used to estimate 
the noise-free versions of i^ and fy as WT^c* and WTyC*, 
respectively. These estimates can be subsequently passed on 
to the fitting procedure of Section III to recover the values of 
(j), which, in combination with a known aperture function A, 
provide an estimate of the PSF i as an inverse discrete Fourier 
transform of the autocorrelation of P = AeJ'^. Algorithm 1 
below summarizes our method of estimation of the PSF. 

where A = A, = Xy. In this form, the problem (|26f is identical Algorithm 1: PSF estimation via DCS 

to (23 I and hence it can be solved by a variety of optimization 
algorithms p6), p7). 

The DCS algorithm augments CCS by subjecting the mini- 
mization in (|26| to an additional constraint which stems from 
the fact that |18| 

dx dy dy dx ' 

which is valid for any two times continuously differentiable 
4i{x^ y). Thus, in particular, the constraint implies the existence 
of two partial differences matrices D^ and Dy which obey 

DJy = Dyix. (28) 

Consequently, if T^. and Ty are the matrices satisfying T^c — 

1) Data: h^, by, and A > 

2) Initialization: For a given transform matrix W and 
matrices/operators '^x, ^y, D^, Dy, T^ and Ty, preset 
the procedures of multiplication by A, A^ , B and B^ . 

3) Phase recovery: Starting with an arbitrary c^^^ and 
p'-''-' = 0, iterate ( 32 1 until convergence to result in an 

Cx and TyC = Cy, respectively, then (31 1 can be re-expressed 
in terms of c as 



^y-^x^ — J-^x-^y^ 

Bc = 0, 


where B :— DyTx — DxTy. Thus, DCS solves the constrained 
minimization problem as given by 

c* = argmin |i||Ac' ~h\\l + A||c'||i| , (31) 

s.t. Be' = 

optimal c*. Use the estimated (full) partial derivatives 
WTxC* and WTyC* to recover the values of (p over H.. 
4) PSF estimation: Using a known aperture function A, 
compute the inverse Fourier transform of P = A eJ^ to 
result in a corresponding ASF h. Estimate the PSF i as 

The estimated PSF can be used to recover the original image 
u from V through the process of deconvolution as explained 
in the section that follows. 

V. Deconvolution 

The acquisition model ([TJ can be rewritten in an equivalent 
operator form as given by 

V ~ H{u} + V, 


A solution to ( 3 1 1 can be found, for instance, by means of 
the Bregman algorithm L28J, in which case c* is approximated 
by a stationary point of the sequence of iterations produced 

^(t+i) ^ 

arg miiic' 

(211^'- - 1^112- 



I Be' 




where Ti, denote the operator of convolution with the estimated 
PSF i. Note that, in this case, the noise term v accounts for 
both measurement noise as well as the inaccuracies related to 
estimation error in i. 

The deconvolution problem of finding a useful approxima- 
tion of u given its distorted measurement v can be addressed in 
many way, using a multitude of different techniques pTj-p3). 
In this work, we use the ROF model and recover a regularized 

^In this work, we use 5 = 0.5. 


approximation as a solution of 

u* = argniin <^ -\\n{u} - v\\l + 7 ||u|JTy 

u I Z 


where ||u||tv = J J\Wu\dxdy denotes the total variation 
(TV) semi-norm of u. 

One computationally efficient way to solve ( (34] i is to 
substitute a direct minimization of the cost function in (|34]) 
by recursively minimizing a sequence of its local quadratic 
majorizers |31|. In this case, the optimal solution u* can 
be approximated by the stationary point of a sequence of 
intermediate solutions produced by 





,(t+i) ^ argmin„ {^\\u ~ u.(*)||2 + 7 \\u\\tv} , 


where H* is the adjoint of H and fi is chosen to satisfy /i > 
IIH*?^!!. In this paper, the TV denoising at the second step 
of ([35| has been performed using the fixed-point algorithm 
of Chambolle fTTl. The convergence of ( |35] l can be further 
improved by using the same FISTA algorithm of pT| . The 
resulting procedure is summarized below in Algorithm 2. 

Algorithm 2: TV deconvolution using FISTA 

1) Initialize: Select an initial value m'^^; set y'"^ = u^^^^ 
and t(o) = 1 

2) Repeat until convergence: 



y^*'> + ^lH* {v - H{y^*^}} 

(*+i) = 0.5 (1 + v'1 + 4(t(*))2) 


In summary. Algorithms 1 and 2 represent the essence of 
the proposed algorithm for hybrid deconvolution of optical 
images. The next section provides experimental results which 
further support the value and applicability of the proposed 

VI. Results 

To demonstrate the viability of the proposed approach, its 
performance has been compared against reference methods. 
The first reference method used a dense sampling of the phase 
(as it would have been the case with a regular SHI), thereby 
eliminating the need for a CS-based phase reconstruction. The 
resulting method is referred below to as the dense sampling 
(DS) approach. Second, to assess the importance of incorpora- 
tion of the cross-derivative constraints, we have used both CCS 
and DCS for phase recovery. In what follows, comparative 
results for phase estimation and subsequent deconvolution are 
provided for all the above methods. 

A. Phase recovery 

To assess the performance of the proposed and reference 
methods under controllable conditions, simulation data was 
used. The random nature of atmospheric turbulence necessi- 
tated the use of statistical methods to model its effect on a 




Fig. 3. An example of a simulated phase (j> (subplot (a)) along with its partial 
derivatives w.r.t. x (subplot (b)) and y (subplot (c)). 

wavefront propagation. Specifically, in this paper, the effect of 
atmospheric turbulence has been described using the modified 
Von Karman PSF model |34|. A typical example of a GPF 
phase (p is shown in subplot (a) of Fig. [5] In the shown case, 
the size of the phase screen was set to be equal to 10 x 10 cm, 
while the sampling was performed over a 128 x 128 uniform 
lattice (which would have corresponded to the use of 16384 
lenses of a SHI). The partial derivatives d(j)/dx and d(j)/dy 
are shown in subplots (b) and (c) of Fig. [3] respectively. 

In the present paper, the subsampling matrices '^^ and ^y 
were obtained from an identity matrix / through a random 
subsampling of its rows to result in a required compression 
ratio r (to be specified below). To sparsely represent the partial 
derivatives of (j)^ ^ was defined to correspond to a four- 
level orthogonal wavelet transform using the nearly symmetric 
wavelets of Daubechies with five vanishing moments |35J. 

To demonstrate the value of using the cross derivative con- 
straint for phase reconstruction, the CCS and DCS algorithms 
have been compared in terms of the mean square errors (MSE) 
of their corresponding phase estimates. The results of this 
comparison are summarized in Fig. |4]for different compression 
ratios (or, equivalently, (sub)sampling densities) and SNR = 40 

As expected, one can see that DCS results in lower values 
of MSE as compared to CCS, which implies a higher accuracy 
of phase reconstruction. Moreover, the difference in the per- 
formances of CCS and DCS appears to be more significant for 
lower sampling rates, while both algorithms tend to perform 
similarly when the sampling density approaches the DS case. 
Specifically, when the sampling density is equal to r = 0.3, 
DCS results in a ten times smaller MSE as compared to 
the case of CCS, whereas both algorithms have comparable 
performance for r — 0.83. This result characterizes DCS as 


Q H n 

— [] 

— ei 

Fig. 4. MSB of phase reconstruction obtained with different methods as a 
function of r. Here, the dashed and solid lines correspond to CCS and DCS, 
respectively, and SNR is equal to 40 dB. 

a better performer than CCS in the case of relatively low 
(sub)sampling rates. 

Some typical phase reconstruction results are shown in 
Fig.|5] whose left and right subplots depict the phase estimates 
obtained using the CCS and DCS algorithms, respectively, for 
the case of r = 0.5. To visualize the differences more clearly, 
error maps for both methods are shown in subplot (c) and 
(d)(note for better visualization error map values are multiplied 
by factor of 10.). A close comparison with the original phase 
(as shown in subplot (a) of Fig. [3]) reveals that DCS provides 
a more accurate recovery of the original (f), which further 
supports the value of using the cross-derivative constraints. In 
fact, exploiting these constraints effectively amounts to using 
additional measurements, which are ignored in the case of 

To investigate the robustness of the compared algorithms 
towards the influence of additive noises, their performances 
have been compared for a range of SNR values. The results 
of this comparison are summarized in Fig. [6] Since the cross- 
derivative constraints exploited by DCS effectively restrict 
the feasibility region for an optimal solution, the algorithm 
exhibits a substantially better robustness to the additive noise 
as compared to the case of CCS. This fact represents another 
beneficial outcome of incorporating the cross-derivative con- 
straints in the process of phase recovery. 

It should be taken into account that, although the shape of 
does not change the energy of the PSF i, it plays a crucial role 
in determining its spatial behaviour In the section that follows, 
it will be shown that even small inaccuracies in reconstruction 
of (j) could be translated into dramatic difference in the quality 
of image deconvolution. 

B. Image deconvolution 

As a next step, the phase estimates obtained using the CCS- 
and DCS-based methods for r — 0.5 were combined with the 
aperture function A to result in their respective estimates of the 
PSF i. These estimates were subsequently used to deconvolve 



Fig. 5. (Subplot (a)) Phase reconstructed obtained by means of CCS for 
SNR = 40 dB and r = 0.5; (Subplot (b)) Phase reconstructed obtained by 
means of DCS for the same values of SNR and r.; (Subplot (c) and (d)) 
Corresponding error maps for CCS and DCS. 

Fig. 6. MSE of phase reconstruction obtained with different methods as a 
function of SNR. Here, the dashed and solid lines con'espond to CCS and 
DCS, respectively, and r = 0.5. 

a number of test images such as "Satellite", "Saturn", "Moon" 
and "Galaxy". All the test images were blurred with an original 
PSF, followed by their contamination with additive Gaussian 
noise of different levels. As an example. Fig. p\ shows the 
"Satellite" image (subplot (a)) along with its blurred and noisy 
version (subplot (b)). 

Using the PSF estimates, the deconvolution was carried out 
using the method detailed in 1 1 1 1. For the sake of comparison, 
the deconvolution was also performed using the PSF recovered 
from dense sampling (DS) of (p. Note that this reconstruction 
is expected to have the best accuracy, since it neither involves 
undersampling nor requires a CS-based phase estimation. 
All the deconvolved images have been compared with their 




Fig. 7. Satellite image (subplot (a)) and its blurred and noisy version (subplot 



Fig. 8. (Subplot (a)) Image estimate obtained with the CCS-based method 
for phase recovery (SSIM = 0.917); (Subplot (b)) Image estimate obtained 
with the DCS-based method for phase recovery (SSIM = 0.781). 

has been shown to yield phase estimates of substantially better 
quality as compared to the case of CCS. 

In this paper, our main focus has been on simplifying 
the structure of the SHI through reducing the number of 
its wavefront lenses, while compensating for the effect of 
undersampling by using the theory of CS augmented by the 
cross-derivative constraint. The solution was computed using 
the Bregman algorithm, which provides a computationally ef- 
ficient framework to carry out the constrained phase recovery. 
Moreover, the resulting phase estimates were used to recover 
their associated PSF, which was subsequently used for image 
deconvolution. It was shown that the DCS-based estimation of 
(j) with r = 0.3 results in image reconstructions of the quality 
comparable to that of DS, while substantially outperforming 
the results obtained with CCS. 

While the proposed method offers a practical solution to 
the problem of phase estimation in adaptive optics, some 
interesting questions about the theoretical aspects of DCS still 
lay open. In particular, the question of theoretical performance 
of CS in the presence of side information on the source signal 
needs to be addressed through future research. 


This work was supported in part by the Natural Sciences 
and Engineering Research Council of Canada and in part by 
Ontario Early Researcher Award program, which are gratefully 
acknowledged. The authors would also like to acknowledge 
Sudipto Dolui for his helpful comments as well as for provid- 
ing deconvolution codes. 

original counterparts in terms of PSNR as well as of the 
structural similarity index (SSIM) of |36|, , which is believed 
to be a better indicator of perceptual image quality |37|. The 
resulting values of the comparison metrics are summarized in 
Table 1, while Fig.lSlshows the deconvolution results produced 
by the CCS- and DCS-based methods. 

The above results demonstrate the importance of accurate 
phase recovery, where even a relatively small phase error can 
have a dramatic effect on the quality of image deconvolution. 
Under such conditions, the proposed method produces image 
reconstructions of a superior quality as compared to the case of 
CCS. Moreover, comparing the results of Table 1, one can see 
that DS only slightly outperforms DCS in terms of PSNR and 
SSIM, while in many practical cases, the difference between 
the performances of these methods are hard to detect visually. 

VII. Discussion and conclusions 

In the present paper, the applicability of DCS to the problem 
of reconstruction of optical images has been demonstrated. It 
was shown that, in the presence of atmospheric turbulence, 
the phase <j) of the GPF P = AeJ'^ is a random function, 
which needs to be measured using adaptive optics. To simplify 
the complexity of the latter, a CS-based approach has been 
proposed. As opposed to CCS, however, the proposed method 
performs phase reconstruction subject to an additional con- 
straint, which stems from the property of V0 to be a potential 
field. The resulting algorithm (referred to as the DCS method) 


[1] J. Yang, J. Wright, T. S. Huang, and Y. Ma, "linage super-resolution 
via sparse representation," IEEE Transactions on Image Processing, vol. 
19, pp. 2861-2873, 2010. 

[2] M. Elad and M. Aharon, "Image denoising via sparse and redundant 
representations over learned dictionaries," IEEE Transactions on Image 
Processing, vol. 17. pp. 3736-3745, 2006. 

[3] IJ. Mairal, G. Sapiro, and M. Elad, "Learning multiscale sparse 
representations for image and video restoration," Multiscale Modeling 
and Simulation, vol. 7, pp. 214-241, 2008. 

[4] R. T. Paul, "Review of robust video watermarking techniques," IJCA 
Special Issue on Computational Science, , no. 3, pp. 90-95, 2011. 

[5] G. D. Boreman, Modulation Transfer Function in Optical and Electro- 
Optical Systems, SPIE Optical Engineering Press, Bellingham, Wash- 
ington, 2001. 

[6] R. T. Paul, "Blind deconvolution via cumulant extrema," IEEE Signal 
Processing Magazine, , no. 3, pp. 2A-A2, 1996. 

[7] D. Kundur and D. Hatzinakos, "Blind image deconvolution," IEEE 
Signal Processing Magazine, , no. 3, pp. 43-64, 1996. 

[8] S. Geman and D. Geman, "Stochastic relaxation, gibbs distribution and 
the bayesian restoration of images," IEEE Trans. Pattern Analysis and 
Machine Intelligence, vol. PAMI-6, pp. 721-741, 1984. 

[9] V. Ton'e T. Poggio and C. Koch, "Computational vision and regulariza- 
tion theory," Nature, vol. 317, pp. 314-319, 1985. 
[10] L. I. Rudin, S. Osher, and E. Fatemi, "Nonlinear total variation based 
noise removal algorithms," Phys. D, vol. 60, pp. 259-268, November 
[11] A. Chambolle, "An algorithm for total variation minimization and 
applications," Journal of Mathematical Imaging and Vision, vol. 20, 
pp. 89-97, 2004. 
[12] O. Michailovich and A. Tannenbaum, "Blind deconvolution of medical 
ultrasound images: Parametric inverse filtering approach," IEEE Trans- 
actions on Image Processing, vol. 16, no. 12, pp. 3005-3019, December 
[13] W. H. Richardson, "Bayesian-based iterative method of image restora- 
tion," / Opt. Soc. Am. A, vol. 62, no. 1, pp. 55-59, 1972. 








Noise std 








0.005 1 10-^ 








PSNR comparison (in 






































































SSIM comparison 





































































[14] L. B. Lucy, "An iterative technique for the rectification of observed 
distributions," Astron. /, vol. 79, no. 6, pp. 745-754, 1974. 

[15] P. A. Jansson, Deconvolution of Images and Spectra, Opt. Eng. 36, 
3224, 1997. 

[16] D. Dayton, B. Pierson, B. Spielbusch, and J. Gonglewski, "Atmospheric 
structure function measurements with a shack-hartmann wave-front 
sensor," Journal of Mathematical Imaging and Vision, vol. 20, pp. 89- 
97, 2004. 

[17] R. G. Lane and M. Tallon, "Wave-front reconstruction using a shack- 
hartmann sensor," Applied Optics, vol. 31, pp. 6902-6908, 1992. 

[18] M. Hosseini and O. Michailovich, "Derivative compressive sampling 
with application to phase unwrapping," in Proceedings of EUSIPCO, 
Glasgow, UK, August 2009. 

[19] M. C. Roggemann and B. M. Welsh, Imaging Through Turbulence, CRC 
Press, Boca Raton, FL, 1996. 

[20] D. L. Fried, "Statistics of a geometric representation of wavefront 
distortion," J. Opt. Soc. Am., vol. 55, pp. 1427-1431, 1965. 

[21] O. Michailovich and A. Tannenbaum, "A fast approximation of smooth 
functions from samples of partial derivatives with application to phase 
unwrapping," Signal Processing, vol. 88, pp. 358-374, 2008. 

[22] Y. Tsaig and D. L. Donoho, "Compressed sensing," IEEE Trans. Inform. 
Theory, vol. 52, pp. 1289-1306, 2006. 

[23] J. Primot, G. Rousset, and J. C. Fontanella, "Deconvolution from wave- 
front sensing: a new technique for compensating turbulence-degraded 
images," / Opt. Soc. Am., vol. 7, pp. 1598-1608, 1990. 

[24] C. E. Shannon, "Communication in the presence of noise," Proceedings 
of the IRE, vol. 37, no. 1, pp. 10-21, 1949. 

[25] E. J. Candes, J. Romberg, and T. Tao, "Robust uncertainty principles: 
exact signal reconstruction from highly incomplete frequency informa- 
tion," IEEE Trans, hifo. Theory; vol. 52, no. 2, pp. 489-509, 2006. 

[26] D. L. Donoho and Y. Tsaig, "Fast solution of ii-norm minimization 
problems when the solution may be sparse," 2006. 

[27] W. Yin, S. Osher, D. Goldfarb, and J. Darbon, "Bregman iterative al- 
gorithms for ^ J -minimization with applications to compressed sensing," 
SIAM J. Imaging Sciences, vol. 1, no. I, pp. 143-168, 2008. 

[28] J. Cai, S. Osher, and Z. Shen, "Split bregman methods and frame based 
image restoration," Multiscale Modeling & Simulation, vol. 8, no. 2, pp. 
337-369, 2009. 

[29] S. S. Chen and D. L. Donoho, "Atomic decomposition by basis pursuit," 
Siam Journal on Scientific Computing, vol. 20, 1998. 

[30] I. Daubechies, M. Defrise, and C. D. Mol, "An iterative thresholding 
algorithm for linear inverse problems with a sparsity constraint," Com- 
munications on Pure and Applied Mathematics, vol. 57, pp. 1413-1457, 

[31] A. Beck and M. Teboulle, "A fast iterative shrinkage-thresholding 
algorithm for linear inverse problems," SIAM Journal on Imaging 
Sciences, vol. 2, pp. 183-202, 2009. 

[32] A. Savitzky and M. J. E. Golay, "Smoothing and differentiation of 
data by simplified least squares procedures," Anal. Chem., vol. 36, pp. 
1627-1639, 1964. 

[33] A. N. Tikhonov and V. Y. Arsenin, "Solutions of ill-posed problems," 

[34] J. D. Schmidt, Numerical .simulation of optical wave propagation with 
examples in MATLAB, SPIE, Washington, 2010. 

[35] I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Reg. Conf. Series 
in Applied Math. SIAM, 1992. 

[36] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. R Simoncelli, "Image 
quality assessment: From error visibility to structural similarity," IEEE 
Trans. Image Processing, vol. 13, no. 4, pp. 600-612, 2004. 

[37] Z. Wang and A. C. Bovik, "Mean squared error: love it or leave it? - a 
new look at signal fidelity measures," IEEE Signal Processing Magazine, 
vol. 26, no. 1, pp. 98-117, 2009.